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In  this  report  the  organization  and  operation  of  the  FORTRAN  computer 
code  lONSCNT,  which  calculates  average  ionospheric  amplit,  Je  and  phase 
scintillation  conditions  on  a global  basis,  is  described  together  with 
the  auxiliary  statistics  program  DIST.  lONSCNT  is  a completely  revised 
version  of  an  earlier  computer  code  (Fremouw  and  Rino,  1975,  1976).*  A 
detailed  discussion  of  the  theory  and  experimental  data  base  that  was  used 
in  developing  the  new  code  is  presented  in  Quarterly  Technical  Report  7 
on  this  contract  (Fremouw  et  al.,  1977). 

As  with  the  earlier  version  of  the  code,  lONSCNT  is  built  around  an 
empirically  derived  global  rms  electron  density  model.  The  user  specifies 
the  signal  source  and  receiver  locations  together  with  the  day  number, 
time,  sunspot  number,  and  planetary  magnetic  index  K . Any  input  param- 
eter can  be  varied  with  a predetermined  Increment  and  range.  In  addition, 
to  simulate  an  orbiting  satellite  the  user  of  lONSCNT  can  vary  the  source 
location  along  a nominal  circular  orbit  between  two  specified  points  on 
the  earth.  This  provides  a reasonable  approximation  to  any  satellite 
orbit  over  a short  distance. 

The  principal  output  of  lONSCNT  consists  of  a tabular  list  of  phase 
and  amplitude  scintillation  indices  together  with  a set  of  parameters 
that  can  be  used  in  OIST  to  compute  the  complete  first-order  signal 
statistics.  In  addition,  lONSCNT  will  generate  intensity  and  phase  auto- 
correlation functions  for  a selected  subset  of  the  data  points  computed 
during  any  run.  Finally,  lONSCNT  will  calculate  the  Intensity  correlation 
function  for  signals  at  two  separated  frequencies. 

In  the  remainder  of  this  report  we  shall  present  a detailed  descrip- 
tion of  lONSCNT  and  Instructions  for  its  operation  together  with  examples 
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References  are  listed  at  the  end  of  this  report. 
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to  illustrate  the  various  program  options.  A user  interested  only  in 
the  program  operation  may  wish  to  skip  to  Section  III.  For  convenience 
the  mathematical  computations  that  are  performed  in  each  subroutine  are 
described  in  the  Appendix.  The  auxiliary  program  DIST  is  described  in 
Section  IV. 
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II  DESCRIPTION  OF  lONSCNT 


A.  General 

A general  flow  diagram  for  lONSCNT  is  shown  in  Figure  1.  From  the 
specified  locations  of  the  signal  source  and  receiver,  the  location  of 
the  ionospheric  penetration  point  at  the  height  of  the  irregularity  layer 
is  computed  together  with  the  propagation  angles  relative  to  the  local 
geomagnetic  coordinate  system  shown  in  Figure  2. 

The  rms  electron  density  model  effectively  gives  the  intensity  of 
the  ionospheric  turbulence  for  the  particular  location,  season,  time,  and 
magnetic  activity  level.  The  rms  electron  density  model  also  determines 
the  values  of  the  anisotropy  parameters,  a,  b,  6,  and  the  outer  scale 
size  for  the  focus  component,  a^.  The  outer  scale  for  the  scatter  com- 
ponent, a^,  is  scaled  from  (see  Appendix,  Section  1).  The  separation 

scale  8 , is  a fixed  fraction  of  a.  such  that  8,  > a . 
t f f s 

Separate  diffraction  calculations  are  performed  for  the  focus  com- 
ponent in  which  a near-zone  approximation  is  made  and  for  the  scatter 
component  in  which  a full  diffraction  calculation  is  performed.  The  full 
diffraction  calculation  is  very  inefficient  for  large  scale  sizes.  Thus, 
a near-zone  approximation  is  used  for  the  focus  component  to  speed  up 
the  calculation. 

The  outputs  of  the  separate  scatter  and  focus  diffraction  calcula- 
tions are  combined  to  generate  the  first-order  moments  (I,  in  Figure  1) 
for  all  data  points.  For  a selected  subset  of  the  data  points,  the 
second-order  moments  (II,  in  Figure  2)  for  phase  and  intensity  can  be 
computed.  Two-frequency  correlation  functions  for  intensity  are  com- 
puted with  special  calls  to  the  first-order  time  statistics  subroutines. 

To  calculate  the  second-order  autocorrelation  functions,  it  is 
necessary  to  specify  the  drift  rate  of  the  ionosphere  as  well  as  any 
relative  source  motion.  An  ionospheric  drift  rate  value  is  provided  by 
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FIGURE  1 FLOW  DIAGRAM  OF  lONSCNT 


the  global  ms  electron  density  model.  The  user  can  modify  the  model 
value  If  he  so  desires.  Representative  satellite  velocities  are  com- 
puted whenever  the  nominal  orbit  code  Is  utilized. 

The  most  time-consuming  operation  Is  the  scatter-component  diffrac- 
tion calculation.  A full  two-dimensional  numerical  integration  Is  per- 
fomed  for  each  data  point.  The  actual  Integration  Is  done  In  polar 
coordinates  with  the  integration  over  the  angular  coordinate  perfomed 
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only  once  for  a fixed  set  of  values  over  the  radial  coordinate.  There- 
after, spline  Interpolation  Is  used  to  generate  the  Integrand  for  the 
Integration  over  the  radial  variable  (see  Appendix  A,  Section  2). 

The  same  Integration  Is  performed  for  both  the  first-  and  second- 
order  statistics.  For  the  latter,  however,  the  separation  variable  Is 
nonzero,  and  the  Integration  converges  much  more  slowly.  Autocorrelation 
function  values  that  have  not  converged  after  a specified  number  of  terms 
are  marked  with  an  asterisk  on  the  output  list. 


AH  the  subroutines  used  In  lONSCNT  are  listed  in  Table  1 together 
with  a brief  description  of  their  function.  The  general  program  flow  is 
shown  in  diagram  form  in  Figure  3.  All  of  the  output  for  each  data  set 
'j  is  listed  at  the  end  of  the  computation,  where  intermediate  results  are 

• stored.  The  computations  for  the  first-order  moments,  the  second-order 

j;  moments,  and  the  two-frequency  correlation  function  are  executed  under 

:<■ 

I control  of  separate  subroutines  that  set  up  the  necessary  inputs  for  the 

calls  to  the  common  subroutines. 

A flow  diagram  of  lONSCNT  by  subroutines  is  shown  in  Figure  A. 

The  diagram  together  with  the  subroutine  descriptions  in  Table  1 are 
self  explanatory.  The  detailed  mathematical  manipulations  that  are 
involved  are  listed  in  the  Appendix.  The  program  Itself  is  written  in 
' * standard  FORTRAN  IV  language,  and  it  should  run  efficiently  on  any 

j medium-  to  large-scale  computer. 
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Table  1 


SUBROUTINES  USED  IN  lONSCNT 


Name 


Def Initlon 


"'NSCNT 

READIN,  with 
ENTRY  INC 

FINDORB 

SRCHT 

SRCHH 

SLAT 

SLON 

VXYZ 

ANCSD 

MULPRM 

SCINTl 

COORD 

AZNGCA 

MGFLD 


RMSDN 


KINDZ 


SIGCOMl 


SIGF 


RFOC 


Executive  routine. 

Performs  program  Input  and  increments  the 
changing  parameter (s) 

Finds  circular  orbit  at  a given  height  through 
two  points. 

Finds  orbital  phase. 

Finds  hour  angle  of  orbital  plane. 

Calculates  satellite  latitude. 

Calculates  satellite  longitude. 

Calculates  scan  velocity. 

Plane  trigonometry  routine. 

Sets  and  calculates  model  parameters. 

Sets  up  inputs  and  calls  to  calculate  first- 
order  signal  nK>ments. 

Finds  point  2 given  point  1 and  azlmutli  and 
great-circle  angle. 

Finds  azimuth  and  great-circle  angle  from 
point  1 to  point  2. 

Calculates  the  geomagnetic  coordinates  of,  and 
the  field  components  at,  a specified  point  using 
an  earth-centered,  axially  tipped  dipole  mmlel 
of  the  geomagnetic  field. 

Contains  the  worldwide  model  for  strength  of 
F-layer  irregularities  used  in  the  scintillation 
calculation. 

Computes  the  reduced  distance  from  scattering 
layer  to  receiver  for  use  in  evaluating  Fresnol- 
zone  radius. 

Provides  calls  for  computation  of  the  "six  sigmas," 
S,,  and  PRM.S. 

Computes  complex  field  moments  for  focus 
component . 

Computes  higher-order  derivatives  of  phase 
autocorrelation  function  for  focus  component. 
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Name 


SIGS 


ABC 

BRBIVL 


C0EFSP2 

F 


CP  RMS 


PHAS,  with 
ENTRY  PRT 


ERF 


FI,  with 
ENTRY  F2 


TWOFREQ 


SCINT2 


ELPSPRM 

ELPRJCT 

SIGCOM2 

SIMPA 


SPLINE 

BESR 


AS  IN 
ACOS 
ORDRIOP 
ORI)R20P 


Table  1 (Concluded) 


Definition 


Computes  complex  field  moments  for  scatter 
component . 


Computes  anisotropy  coefficients  for  phase  ACF. 

Integrates  asymmetry  correlation  function  for 
scatter  component. 

Computes  spline  coefficients  for  G(u). 


Computes  Integrand  for  polar  integration  over 
ray  variable  and  angle  variable. 


Computes  integrand  for  polar  integration  over 
ray  variable  from  spline  coefficients. 


Controls  calculation  of  rms  phase. 


Evaluates  the  Hatfield  phase  distribution  of  a 
scintillating  signal. 


Computes  the  error  function. 


Provides  Integral  functions  for  evaluating 
moments  of  phase  distribution. 


Sets  up  calling  sequence  for  two-frequency 
correlation  computation. 


Sets  up  inputs  and  calculates  second-order 
moments . 


Calculates  axial  ratio  and  orientation  angle  for 
spatial  autocorrelation  function. 


Transforms  coordinates  from  plane  normal  to 
propagation  vector  to  receiver  plane. 


Provides  calls  for  computation  of  second-order 
moments . 


Adaptively  performs  numerical  integration  using 
Simpson's  rule. 


Calculates  coefficients  for  spline  interpolation. 


Computes  the  K Bessel  Function  for  a given 
argument  and  order. 


Inverse  sine  routine. 

Inverse  cosine  routine. 
Outputs  first-order  results. 
Outputs  second-order  results. 
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FIGURE  4 lONSCNT  FLOW  DIAGRAM  BY  SUBROUTINES 
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Ill  INSTRUCTIONS  FOR  USE  OF  lONSCNT 

A.  Program  Inputs 

To  run  lONSCNT  it  is  necessary  to  prepare  three  types  of  input 
cards;  (1)  a title  card,  (2)  namelist  cards,  and  (3)  a changing- 
parameter  (s)  /options  card.  Each  type  of  card  is  described  below.  Three 
sample  sets  of  input  cards  are  shown  in  Figure  5. 

The  title  card  is  used  to  assign  any  desired  title  of  up  to  80 
characters  to  a specific  case  or  series  of  cases.  It  is  read  by  the 
executive  program,  lONSCNT,  and  printed  on  the  output.  It  has  no  effect 
on  the  calculation.  Any  number  of  cases  can  be  run  under  a common  title. 

The  namelist  cards  are  used  to  input  the  fixed  parameters  for  a 
calculation  case.  They  are  read  by  subroutine  READIN,  by  means  of  the 
namelist  INPAR.  The  user  gives  the  name  and  desired  value  for  each 
quantity  listed  above  the  dashed  line  in  Table  2 that  is  to  remain  fixed. 
The  number  under  the  heading  K in  the  table  is  a program  index  used  to 
identify  the  parameters.  The  user  may  list  the  fixed  parameters  in  any 
order  on  the  namelist  cards.  A parameter  input  on  a namelist  card  will 
keep  the  specified  value  until  changed  by  means  of  a new  namelist  card 
or  until  input  on  a changing  parameter (s) /options  card. 

The  first  column  of  each  namelist  card  must  be  blank.  Columns  2 
through  7 of  the  first  entry  in  a sequence  must  be  $INPAR.  All  namelist 
cards  except  the  last  must  end  with  a constant  followed  by  a comma.  The 
last  in  a sequence  must  terminate  with  $. 

The  changing-parameter (s) / options  card  contains  three  types  of 
Information:  (1)  the  Independent  variable (s)  for  a given  calculation 

case,  (2)  the  second  frequency  for  the  two-frequency  correlation  calcula- 
tion, and  (3)  the  information  for  second-order-statlstlcs  calculations. 

If  neither  (2)  nor  (3)  is  desired,  this  card  returns  to  the  form  of  the 
changing-parameter (s)  card  of  previous  versions  of  the  program. 
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Table  2 

DATA-CARD  VARIABLES 


K. 

Variable  Name 

Definition 

■■■ 

FREQ 

Frequency  (MHz) 

SSN 

Sunspot  number 

DAY 

Day  of  year 

TIME* 

Meridian  time  at  the  receiver  (decimal  hours) 

RLAT 

Receiver  latitude  (decimal  degrees) 

RLON 

Receiver  longitude  (decimal  degrees)  (+east) 

HR 

Receiver  altitude  (meters) 

8 

TLAT 

Transmitter  latitude  (decimal  degrees) 

■ 

TLON 

Transmitter  longitude  (demical  degrees) 

(+east) 

HT 

Transmitter  altitude  (meters) 

H 

FRP 

Planetary  magnetic  index,  K (0.0  through 

9.0)  P 

— 

Not  implemented 

H 

11 

Data  flag  ••  0.1,  or  2 (see  text) 

13 

RCRD 

Receiver  coordinates  (decimal  degrees) 

14 

TCRD 

Transmitter  coordinates  (decimal  degrees) 

15 

ORBT 

Transmitter  coordinates  (decimal  degrees) 

Note  that  TIME  has  assumed  a different  definition  from  the  definition 
that  was  used  in  previous  versions  of  the  program.  It  is  now  meridian 
time  at  the  receiver.  In  the  case  where  ORBT  is  the  changing- 
parameter  Indicator,  TIME  is  meridian  time  at  the  receiver  for  the 
first  orbit  point. 
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(See  Table  2 regarding  TIME  exception.)  The  absence  of  additional 
entries  indicates  that  the  added  calculations  are  to  be  suppressed. 

If  the  variable  named  on  the  changing-parameter (s)  card  corresponds 
to  a K value  in  Table  2 between  1 and  12,  the  calculation  will  be  per- 
formed as  a function  of  that  single  changing  parameter.  If  RCRD  (or 
TCRD)  is  specified,  receiver  (or  transmitter)  position  will  be  incre- 
mented along  a great-circle  arc;  thus  the  calculations  will  be  dependent 
upon  the  variation  of  both  a latitude  and  a longitude.  If  ORBT  is 
specified,  the  transmitter  will  move  (latitude  and  longitude)  along  a 
circular  orbit  at  a fixed  altitude  (HT) . 

The  parameters  that  have  been  added  to  the  changing-parameter (s) / 

I options  card  relate  to  the  two-frequency  correlation  calculation  and  to 

1 the  second-order-statistics  calculations.  Their  formats  can  be  seen  in 

\ ' Table  3.  If  FREQ2  » 0.0,  the  two-frequency  correlation  option  is  by- 

i ! passed;  otherwise,  the  value  is  the  second  frequency. 

' I If  NSCND  ^ 0,  the  second-order  calculations  are  bypassed.  NSCND  - 1 

• I gives  second-order  output  for  all  calculation  points;  NSCND  ~ 2 gives  it 

I : for  points  1,  3,  5,  ...;  NSCND  “ 3 gives  it  for  points  1,  4,  7,  ...;  etc. 

Second-order  calculations  are  done  for  11  points  from  time  lag  • 0.0  to 
j DLTMMX  seconds.  If  DLTMMX  ^ 0.0,  only  the  second-order  heading  informa- 

i tion  is  given,  but  no  detailed  second-order  moment  calculations  are  per- 

formed. If  the  changing-parameter  option  is  ORBT,  the  scan  velocity 
! components  of  the  receiver  as  seen  from  the  penetration  point  are  com- 

i puted  and  VXI,  VYI,  and  VZI  (drift  velocity  component)  are  added  to 

I obtain  VX,  VY,  and  VZ  used  in  the  second-order  calculations.  If  this 

I option  is  other  than  ORBT,  the  user  specifies  VX  - VXI,  VY  ■ VYI,  and 

IVZ  ■ VZI.  If  VXI  - VYI  - VZI  ■ 0.0,  the  program  supplies  nominal  values. 

After  a set  of  cards  is  read  and  the  requested  options  are  executed 
and  the  outputs  listed,  the  present  value  of  11  is  examined.  If  this 
[ value  is  0,  the  program  attempts  to  read  a new  title  card;  if  it  is  1, 

the  program  attempts  to  read  a new  namelist  card;  and  if  it  is  2,  the 

! program  terminates. 

( 
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Table  3 ! 


CONTENTS  OF  CHANGING-PARAMETER(S) /OPTIONS  CARD 


Column 

Numbers 

Name  of 
Variable 

Contents 

For  One  Variable 

For  Two  Variables 

1-4 

NAME 

Changing-parameter 

Changing-parameter 

Indicator 

Indicators 

9-16 

Al 

Initial  value 

Initial  latitude 

17-24 

A2 

Final  value 

Initial  longitude 

25-32 

A3 

Increment 

Final  latitude 

33-40 

A4 

— 

Final  longitude 

41-42 

NA5 

— 

Number  of  increments 

Definition 

43-44 

NSCND 

Specifies  calculation  points  on  which  to 

perform  second-order 

calculations. 

45-48 

DLTMMX 

Maximum  time  lag,  In 

seconds . 

49-56 

VXI 

In  ORBT  mode,  velocity  components  added  to 

57-64 

VYT 

components  due  to  source  motion;  otherwise. 

total  velocity  component  used  In  second- 

65-72 

XZI 

order  calculation.  Program  supplies  values 

If  VXI  - VYI  - VZI  - 

0.0  (see  text). 

73-80 

FREQ2 

Second  frequency  (MHz  for  2-frequency 

correlation) . 

1 

For  model  development  and  testing,  a plot  routine  was  used  In  con- 
junction with  lONSCNT  to  produce  graphs  such  as  those  In  Figures  5 
through  11  of  Quarterly  Technical  Report  7 (Fremouw  et  al. , 1977).  The 
available  version  of  lONSCNT  provides  arrays  that  the  user  can  employ 
with  his  own  computer  peripherals.  In  the  absence  of  user  modification, 
the  arrays  will  be  output  In  tabular  form  along  with  Identifying  Informa- 
tion, as  shown  In  the  sample  calculations  below. 

B.  Sample  Calculations 

In  Figure  6 we  show  the  standard  lONSCNT  output  for  a case  In  which 
the  orbit  option  was  used  and  two-frequency  Intensity  correlations  were 
requested.  The  card  Inputs  are  shown  in  Figure  5<a).  The  title  is 
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FIGURE  6 SAMPLE  OUTPUT  FROM  lONSCNT  SHOWING  TWO-FREOUENCY  INTENSITY  CORRELATION 


supplied  by  the  user.  The  fixed  parameters  for  the  run  are  listed  to 
verify  the  input  conditions.  The  six  sigmas  and  the  intensity  and  phase 
scintillation  indices  are  listed  for  each  point  together  with  the  Briggs- 
Parkin  angle  (the  angle  between  the  local  magnetic  field  direction  and 
the  propagation  vector)  and  the  variable  parameters  for  the  run.  The 
points  are  numbered  for  reference. 

In  the  second  list  of  numbers  the  first-order  moment  for  the  second 
frequency  is  listed  together  with  the  two-frequency  intensity  correlation 
for  intensity,  CIF1F2.  The  two-frequency  intensity  correlation  is 
described  in  the  Appendix,  Section  6. 

In  Figure  7(a)  we  show  a second  example  using  the  orbit  option.  For 
this  example,  however,  the  second-order  moments  were  requested  for  every 
seventh  point.  The  data  cards  are  shown  in  Figure  5(b).  The  second- 
order  output  is  shown  in  Figure  7(b). 

The  "effective"  velocity  VEFF  is  defined  in  the  Appendix,  Section  2. 
The  effective  velocity  is  a function  of  the  true  scan  velocity  and  the 
anisotropy  parameters.  It  accounts  for  the  fact  that  for  the  same  scan 
velocity,  the  time  structure  of  the  actual  signal  is  more  rapid  when  the 
scan  is  across  the  principal  irregularity  axis  than  when  it  is  along  the 
principal  irregularity  axis.  The  azimuth  of  the  scan  velocity  is  the 
true  azimuth  of  the  scan  velocity  component  perpendicular  to  the  line  of 
sight  projected  onto  the  receiver  plane  (tangent  to  the  earth). 

The  axial  ratio  and  true  azimuth  of  the  contours  of  constant  cor- 
relation in  the  receiver  plane  are  also  listed.  The  lag  parameter  in 
meters  is  measured  along  the  scan  direction.  The  listed  parameters  are 
shown  schematically  in  Figure  8 for  point  8.  By  comparisi.  i with  the 
corresponding  parameters  for  point  1,  it  can  be  seen  that  for  point  1 the 
scan  velocity  is  closer  to  the  direction  of  the  semi-minor  axis  and  the 
ellipse  is  more  elongated.  Thus,  the  effective  scan  velocity  is  larger. 

The  correlation  functions  for  phase  and  amplitude  are  listed  under 
the  columns  marked  COMPOSITE  [see,  for  example.  Figures  7(b)  and  9(b)]. 
For  reference  purposes  the  separate  phase  and  amplitude  correlation 


19 


mis  fAQ®  IS  BEST  lUAIilTY  PRACJICAW* 

ijBmcoriii'.uiwuisiuaDroj)D,Q  — 


I ' ♦ I 

. »«*.»«  n « 

V n r !•  .1  *•*  iV  iM  ► 4 

« • »•  a t r • i»  •>» 

a « «(  W'  • c u « «.'  aU 
♦ »w»r  wi  «»••>«< 

/••'*•  *1  •••'•  • 


'*•  0|«  .s  • r(M^o«^K*v 

/ #.«  O 4 a^  % r ■!  rv  <4  • • *v  <V 

r-  ,fA  «•  o <1  ««  (V  4*  f A »•  4 

av  C O 1*.  a I * 4 iM  ila  tfi 

»-■  a a.'  '4<  ^ Mi  9 <*l  IW  iW  '■* 

•^S  • •••,•••>•••  • 

- I 


4»  f r i»  4 '•'a^  •'a 

0 Q n n 0 o t\ 

■0  f.  *>  4>  4'  ••>  «n  f\4  4>4  fV 


o -«  o u iT  ■>>  r» 
a3  ji  ri  a4  41  /a  m 'M  <4  r 
0<^*>n4/a«^/l4t 
O K.*4  U>II4  4>u>4.4 
0 ^ K tia  t >*1  IV 


t»  0 ty0%0>n(S40ti^0n 

04>«w-t04<Mafl44»4IO 
Oi>4^  J M»^aN.  V4>r«  W 

a n‘na.aW,/aA4> 

•A  ^ -r  ^ a %,  c 


^^44  4l4i««»^04>0 

«?  t«  n 4t  4«  <*>  o « .0  r\«  I 


•»  a ^ to  0 iT  a^ 

O n 0«n^»-a4^f«)|«^ri 

a«  4 -«  cp  .t1  « 4>  n ap.in  0 

«a  «n  ^ •»  a^  4 • IV  »** 

4 af  « 4»  I'  At  4 mj  •*«  IV 

(•  »•••%•*•{•• 


14  (If  ••! 

J A t»  '*t  0 

J I J a*  > •» 

ff  1)  on  V 

•>,  a on>  ► 

o.u.  •,  * • 


•I  U O 4 <«•  aj 

0“  0-  0^  *O>«*V044 

4 a •«><'«  1*1  4 ^ IV  «* 

••44  a>  lU  an  4 n IV  iV  '4 

tU  *,444 


«V  ( 't  If  M 4 aa>  • 
■•  . • . 4 rv.  rv  Ai  a 1 

t4<'U<4>W4A 
I « IW  4'<  • a O aj  a 


'.»  » il  4<  <tl  n 4^  4 V aJ 

*^*'•.#4<}1l•nn 
. • oi  » . 1 r n V O 4>  tf  I 

a ' f >V  4 At  y 4V  • 4 .k' 

\»  0 i>  4\  0 $ 0 fi  At  A\  4 


« ■•!  aJ  aa,  j «|  .y  0 ,|tl 

aa  u 4 a4  a M U ii  At  aa , 

a a#  vuit>ai^  •aa; 
A «a  t (•  IV  t A 4 ft  ‘VI 

4 IV  '‘4  r-a  O a a a.a  > Oi 


O 4 ai 

*4  A)  -V  ‘«- 

A O l\j  n 'V  •>  '-a  |t  t 0 

i(  v)  at  ’0  ■»•  A'  IV  At  --a  r- 

Oi  t3<A0  0\0an»r» 

t a • • a 

r;  -1 


• 4 a^  » >•  r»  ji  *1 

I 4»  O A*  At  V -4  •*>, 

. ai  (a  u r -n  a > 


a-  ' O ap  « 4 

•4  A»  O ^ 

J Jt  O p At  r* 

tV  iJ  at  •»>  'P  4 

.to  0 » tr 

a ‘ -I  • • ' 


« .k  e «n  1*1  -u 

t aai  O > a > 4 n 

at  ^ O > a4  It 

4«  at  iw  •■4  >ti 

^ at  ^ 4<  p rt  IV  • 


a.4».a^A><»P<->t't4 
t V V U P At  .4  r ‘t  al 

» »'I  »»  n ^ ra.  At  •«  / .P 

»-  4 ‘t  At  A jt  .*  * ^1 

• t a**  P « ■»-  P M« 

•,•••••  t ••  1 


44.  . -a  A*  1*1  P « -a  IV  at  P •# 

IaA  aJ<M44r>4'apa^|>At 

X4-  a.  *t<t/t  vnVPHAt> 

0 a-a  ^ .al  Al  •*•  V • n 4>  4A 

a«  A <•  ••  At  A •*»  V >4  r*  a > 


^IPtA  VlV  PlP^•  ^>ra 
• 4>  Al  0 •!  ’V  .V  O lf> 

4 4 n 4 At  «i»|p  a 


O 0-  <0  »a  1*1  Wt  fa  O 

^ ' V It  0 V At  'W  ^4 

1 -4  -4  -•  aj  V VI  n 


' III  at  .»1  *a  j|  .»  A'  a n O- 

I ^ m .»»  .V  4a  4 p , a ga 

At  a » 41  -4  t.  I,  ,0  f 

4 4 V‘  V •>  •*>  4 4 At 


« <p4>.^««IVAllVn 

I 

I 


FIGURE  7 (Concluded) 


NORTH 


FIGURE  8 DIAGRAM  SHOWING  RELATIVE  ORIENTATION  OF  SCAN  VELOCITY 
AND  CONTOUR  OF  CONSTANT  SPATIAL  AUTOCORRELATION 


functions  for  the  scatter  and  focus  components  are  also  listed.  To 
reiterate,  the  lag  parameter  in  meters  is  a measure  of  the  correlation 
distance  along  the  direction  of  scan. 

As  a final  example,  in  Figures  9(a)  and  9(b)  we  show  the  outputs 
for  a run  with  the  source  and  receiver  fixed  and  frequency  varying  from 
50  to  500  MHz  in  steps  of  50  MHz  with  second-order  outputs  requested 
for  every  fourth  point.  The  corresponding  input  cards  are  shown  in 
Figure  5(c) . 

The  user  is  warned  that  the  full  second-order  outputs  require  con- 
siderable computation.  Thus,  they  are  time  consuming  and  therefore 
expensive  to  generate.  The  first  example  presented  in  this  section  used 
143  s of  CDC  6400  computer  time  including  the  generation  of  plots  not 
shown  here.  By  comparison,  the  second  example  used  ~ 400  s of  CDC  6400 
computer  time. 
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calculated  from  a model  developed  RV  STANFORD  RESEARCH  INSTiTUTCt 
MENLO  PARKt  California. 
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IV  THE  AUXILIARY  STATISTICS  PROGRAM  DIST 


A.  General 

Program  DIST  is  to  be  used  in  conjunction  with  Program  lONSCNT  to 
provide  the  user  with  detailed  fading  statistics  for  a signal  that  has 
passed  through  the  ionosphere.  The  necessary  input  parameters  for  DIST 
are  computed  in  Program  lONSCNT  as  described  in  Sections  II  and  III. 
Program  DIST  provides  the  user  with  statistical  Information  about  the 
received  signal  in  the  form  of  probability  density  functions  for  amplitude 
and  phase  (tables  and  plots).  The  user  can  estimate  fading  statistics 
from  the  cumulative  probability  density  functions  for  amplitude.  A 
cumulative  probability  density  function  of  phase  is  also  provided. 

B.  Description  of  the  Program 

Program  DIST  computes  the  density  distribution  of  phase  and  ampli- 
tude for  the  two-component  model  of  scintillation  statistics.  The  two- 
component  model  can  be  characterized  as  the  product  of  two  statistically 
independent  components  having,  respectively,  Gaussian  and  log-normal 
statistics  as  described  in  Quarterly  Technical  Report  7 (Fremouw  et  al., 
1977).  The  joint  probability  density  function  (PDF)  of  amplitude,  a, 
and  phase,  (|),  is  given  by  the  two-dimensional  integral 
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In  Eq.  (1),  P (...)  the  joint  PDF  for  the  real  and 
X y 
s s 


imaginary  parts  of  the  scatter  component,  while  Py  * (...)  is  the  joint 

XfiPf 


PDF  for  the  log-amplitude  and  phase  of  the  focus  component.  Both  PDFs 


are  assumed  to  have  the  general  joint  Gaussian  form.  Thus,  they  are 


2 2 2 2 

specified  in  terms  of  the  six  sigmas,  a , a , a , a o.  . and  a , , 

X ’ y * xy’  Xf  . 4'  X't' 


which  are  variances  and  covariances  of  the  respective  components.  These 


parameters  are  input  to  this  program  and  are  given  by  Program  lONSCNT 


as:  - SGX2;  - SGY2;  o ■ CRXY;  o ^ - SGLA2;  - SGPH2;  o 

X y xy  ’ Xc  X'P 


f 

CRLAPH.  The  mean  values  of  x and  y for  the  scatter  component  and  x snd 

2 2 2 

for  the  focus  component  are,  respectively,  n ■ 1 - o - a and  0; 

2 X X y 

and  in (a)  + a and  0.  These  mean  values  ensure  energy  conservation  for 


the  composite  signal  and  the  scatter  and  focus  components  separately. 


The  user  is  generally  Interested  in  the  statistics  of  amplitude  and 
phase  separately.  Thus,  we  Integrate  Eq.  (1)  over  the  phase  variable 


to  obtain  the  PDF  for  amplitude,  P^(a),  which  can  be  written  as 


- a 


*.n(a)  + - X„  I dX 


) 


(2) 


In  Eq.  (2),  P.  (•)  is  the  amplitude  PDF  for  the  scatter  component  and 

A 

s 


p (•)  is  the  log  amplitude  PDF  for  the  focus  component. 


'f 


A similar  calculation  gives  the  PDF  for  phase: 


00 

■ f (♦»)  (*  ■ 

«/.oo 


d4> 


(3) 


In  Eq.  (3),  P.  (•)  is  the  phase  PDF  for  the  scatter  component — the 


Hatfield  distribution  (Hatfield  and  Rino.  1975 — and  P (•)  is  the  phase 

^f 

PDF  for  the  focus  component. 
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In  addition  to  the  above  three  distributions,  PDFs  and  P^(4i), 

the  corresponding  cumulative  density  functions  (CDF),  are  computed.  For 
the  amplitude,  the  CDF  is  presented  in  dB  for  fade-margin  computations. 

The  output  represents  the  probability  of  fading  below  a specified  ampli- 
tude. The  CDF  of  phase  represents  the  probability  of  a phase  change 
greater  than  a specified  value.  (Both  positive  and  negative  changes  are 
considered.)  The  Inputs  to  DIST  are  generated  by  lONSCNT. 

Table  4 gives  a list  of  the  subroutines  required,  together  with  a 
brief  description  of  the  function  of  each  routine.  The  operation  of 
DIST  is  straightforward.  The  various  functions  are  evaluated  to  generate 
the  necessary  Integrands,  which  are  then  evaluated  by  numerical  integration 

Table  4 

SUBROUTINES  USED  IN  PROGRAM  DIST 


Subroutine 


COMP 2 A 
FEHLIN 


GAUSS 


GAU2ZM 


GRAPH4 


ITERAT 

NUPHAS 

PAAMP 


Function 


Control  program,  reads  input;  calls  subroutines. 
Computes;  lists;  plots  CDF  of  Eq.  (3). 

Computes;  lists;  plots  CDF  of  Eq.  (2). 

Computes  Eq.  (2). 

Computes  error  function  complement. 

Used  to  compute  limits  of  integration 
for  Eq.  (2). 

Computes  one-dimensional  Gaussian  density 
function. 

Computes  bivariate  Gaussian  zero  mean 
density  function. 

System-dependent  plotting  routine. 

User  may  substitute  his  own  or  a dummy  routine. 


Finds  zero  of  a function  by  iteration. 
Computes  p^(')  and  stores  In  array. 
Computes  integrand  of  Eq.  (2). 

Lists  and  plots  Eq.  (2). 


Reference 


PRTHET 


PGTHET 

PAAMP 

PRAPH 
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Table  4 (Concluded) 

Subroutine  Function 

PDFP  Lists  and  plots  Eq.  (3). 

PGTHET  Computes  Integrand  of  Eq.  (3). 

PNEW  Computes  stores  in  array. 

PR  Computes  P^  (•)• 

s 

PRAPH  Computes  Integrand  needed  to  find  P^  (•). 

s 

PRTHET  Computes  p,  (•) 

*s 

QG5  Five-point  Gausslar  quadrature. 

QG6  Six-point  Gaussian  quadrature. 

QG6L2  Six-point  Gaussian  quadrature  in  two  parts 

QG16  Sixteen-point  Gaussian  quadrature. 

QG20  Twenty-point  Gaussian  quadrature. 

SYGN  Return  plus  or  minus  one,  depending  on  sign 

of  argument. 


Instructions  for  Use 


Reference 


NUPHAS 

C0MP2A 

PAAMP 


PGTHET 


NUPAAS 

ITERAT 


1.  Input 

The  input  required  for  this  program  is  one  card  per  data  set. 

The  card  contains  the  six  parameters  of  the  two-component  distribution 
which  are  obtained  from  the  output  of  lONSCNT.  The  parameters  are  SGX2, 
SGY2,  CRXY,  SGLA2,  SGPH2,  and  CRLAPH.  They  must  be  written  with  a decimal 
point  in  ten-column  fields.  A blank  card  indicates  the  end  of  all  data 
sets. 

The  user  must  substitute  his  own  plotting  subroutine  for  GRAPH4 
if  he  desires  plots.  If  listings  only  are  desired  a dummy  subroutine 
GRAPH4(...)  may  be  supplied. 

2.  Output 

The  output  consists  of  tables  and  plots  of  four  distributions: 
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• P(A),  the  amplitude  probability  density  versus  amplitude 
normalized  to  the  undisturbed  amplitude. 

• CPF,  the  cumulative  amplitude  probability  versus 
signal  strength  in  dB  relative  to  the  undisturbed  level. 

• P(PHI)  the  phase  probability  density  versus  phase 
angle  in  radians. 

• Probability  of  phase  change  greater  than  THETA  in  radians. 

3.  Examples 

We  present  two  examples  to  illustrate  the  use  of  DIST.  The 
card  Inputs  for  both  examples  are  shown  in  Figure  10.  In  the  first 
example  (Figure  11),  very  large  phase  changes  occur,  so  the  cumulative 
phase  density  is  not  meaningful  and  therefore  is  not  computed.  In  the 
second  example  (Figure  12),  all  four  distributions  are  shown.  The  data 
output  format  is  self-explanatory.  The  running  time  on  the  CDC  6400 
for  the  two  examples  was  29.2  sec.  This  time  Includes  computation, 
plotting,  and  listing. 
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FIGURE  11  SAMPLE  OUTPUTS  FROM  DIST  FOR  HIGH  LEVEL  OF  ACTIVITY 
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Appendix 

SUMMARY  OF  THEORETICAL  CALCULATIONS  USED  IN  lONSCNT 


! 

? 

1. 


1.  Phase  Autocorrelation  Function 

The  phase  autocorrelation  used  in  lONSCNT  is  described  in  Rino  and 
Fremouw  (1977).  It  has  the  general  form 

R^(Ap)  - AL  R(f(APg)j  (A-1) 

where 

Ap  » Ap  - tan  0 a.  Az  (A-2) 

S ‘Sp 


and 


tc  - K^_j^y2(26/a)/K^_l(2e/a)  . (A-6) 

The  anisotropy  of  the  irregularities  manifests  itself  in  the  A,  B, 

C coefficients.  These  coefficients  are  given  in  terms  of  the  elongation 
parameters  along  and  transverse  to  the  field  (a  and  b,  respectively),  the 
orientation  angle  of  the  transverse  axis,  6,  the  magnetic  dip  angle,  iji, 
and  the  propagation  angles  6 and  V’ , as  follows: 


A = 

f''  " 2 2 

Cii  + ® ^ 

A 

-^<=13 

tan  9 cosv’j 

(A-6a) 

B - 

f"  " 2 

12  ^ ^33  ® ^ 

COS  <f} 

- tan  0 (Cj^3  sln<p 

(A-6b) 

+ 623  cos  <pj 

(A-6b) 

C - 

f"  " 2 2 

1 C22  + ® 

- 2^23 

tan  9 sin  </>j  .• 

(A-6c) 

where 

^11 

2 2 2 2 
= a cos  li)  + sin  ^ (b 

sin^  6 

...  2 .. 

+ cos  6) 

(A-7a) 

C22 

k2  2.^,2. 

■ b cos  6 + sin  6 

(A-7b) 

*^33 

2 2 2 2 
“ a sin  iji  + cos  ij;  (b 

sin^  6 

+ cos^  6) 

(A-7c) 

^12 

2 

■ (b  - 1)  sin  ip  sin  6 cos  6 

(A-7d) 

A 

"l3 

/ 2 u2  4 2 . 2 

- (a  - b sin  o - cos 

6)  sin 

4'  cos  tp 

(A-7e) 

"23 

2 

■ - (b  - 1)  cos  (p  sin  6 

cos  6 

(A-7f) 
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For  the  scatter-component  calculations  we  have  used  the  limiting 
forms  of  Eqs.  (A-5)  and  (A-6)  as  By  using  the  small-argument  forms 

of  (viz,  K^(2x)  ~ j r(Y)  X ^) , we  have 


R,(y) 


(A-8) 


and 


iCg  - a 2^(Y-i)/r(Y-l) 


(A-9) 


Equations  (A-1)  through  (A-7)  together  with  the  limiting  forms  of  Eqs. 
(A-8)  and  (A-9)  constitute  a complete  characterization  of  the  Integrated 
phase  along  the  principal  raypath,  which  Is  the  starting  point  for  phase- 
screen  calculations. 

The  spectral  density  function  for  phase  has  the  form 


2 2 2 
r A AL  sec 
e 


-tan  6 a, 


(A-10) 


where 


-tan  0 a, 

4 


ab<AN^>Q 


+ BA  K 
X y 


(A-11) 
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Q(q)  - (2ii  a0) 


5 -<r>-4)  K 


V4t  ! f^),. 

s-1  !) 


(A-12)' 


and  (KfK^)  is  the  three-dimensional  spectral  density  function  of  the 

0 

electron  density  Irregularities.  When  a >>  3 and  2v/a  < q < 27t/3, 


3/2  r(Y  + 1/2)  3 

Q(q)  « (2ti)  ~r(Y  r i)"  “ 


[1  + (aqr/2] 


2 /oi Y+I/2 


(A-13) 


1 -4 


In  the  model  we  have  taken  y - 1.5,  so  that  Q(q)  « a q . It  fol- 
-1  -4 

lows  that  a q . For  the  focus  component  we  take 


<4»2>  - , 


(A-14) 


which  makes  the  composite  spectral  density  function  essentially  a con- 
tinuous power  law  except  where  2tt/3j  ^ q < 

In  the  computations,  the  temporal  correlation  functions  are  obtained 
by  replacing  Ap  and  Az  in  Eq.  (A-2)  by  the  relative  drift  velocities  v^, 

V , and  V . We  have  also  used  the  definition 

y 2 


^ff  ■ '<"s)  • 


(A-15) 


2.  Diffraction  Calculations  for  the  Scatter  Component 

In  the  Gaussian  phase-screen  model  we  can  easily  calculate  the 

correlation  functions  R (Ap)  - <uu  *>  and  B (Ap)  “ <uu  > at  the  edge  of 

s s 

the  phase  screen.  Here  u is  the  scalar  wave  field.  The  results  are: 


44 


R^(ip) 


■ exp  0^  AL  (1  ± R^(Ap)]|  - exp  Al|  (A-lb) 


where  the  upper  sign  Is  used  for  R and  the  lower  sign  for  B 

s s 


The  mutual  coherence  function  R^  is  propagation-invariant,  whereas 


for  2 > 0, 


B^(Ap;z)  ■//  B^(Ap' ) F(^  - a'p' ) d&p  (A-17a) 


where 


Fit)  » 


cos  0 + (MO‘ 


(A-17b) 


(cos^  sin^\ 
-sin<p  cos  <p  / 


(A-17c) 


By  making  appropriate  substitutions  and  several  variable  changes, 
Eq.  (A-16a)  can  be  written  in  the  form 


|®^s|  2 exp{-  al  Al}/*”  , 2 


® (A-18) 


where 


G(w;yj,) 


(A-19) 


r 


A5 


(A-27) 


tan  2<( 

K 


In  the  numerical  Integration  performed  In  BRBIVL,  G(w;yjj)  Is  first 
evaluated  using  an  adaptive  Simpson  method  for  a logarithmically  spaced 
set  of  u values.  A set  of  spline  coefficients  Is  then  generated  for 
these  points.  For  the  u Integration,  G((i);y^)  Is  Interpolated  using  the 
spline  coefficients. 

The  Integration  over  the  u)  variable  Is  performed  by  Integrating 
over  a number  of  cycles  of  the  sin  or  cos  terms  and  then  summing  the 
resultant  series.  The  Integration  Is  terminated  when  the  resultant 
series  converges  to  within  a preassigned  value. 


3.  Propagation  Calculations  for  the  Focus  Component 

For  the  focus  component,  the  relevant  quantities  are  the  log- 
amplitude  and  phase  of  the  field.  These  are  calculated  by  using  the 
ray-optics  approximation 


aXf 

at 


k ’"i*t 


(A-28) 


from  which  It  follows  that 


<Xf't'f 


a^R, 


a^R, 


a^R 


’ll 


+ 2a 


12  aAp  aAp 

X y 


+ a 


1 


22  0Ap2 

y 


(A-29) 
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and 


<XfXf 


(\^  sec  0\ 

>■ 


*11 


aAp 


r * '‘*11*12 


aAp  8Ap 
X y 


(^*12  ^ ^*11*22) 


8^, 


a^R. 


a^R, 


aAp^8Ap^ 
X y 


^®12®22  2„.  3 ^22  4 


8Ap  8Ap' 
X ^y 


aAp, 


(A-30) 


The  evaluation  of  the  derivatives  Is  tedious  but  straightforward. 

The  results  are  summarized  In  Figure  A-1,  which  Is  a partial  listing  of 
the  subroutine  RFOC.  The  variable  names  In  the  code  are  self  explanatory. 

4.  Calculation  of  Intensity  and  Phase  Moments  for  Composite  Signal 

The  phase  autocorrelation  function  Is  a simple  superposition  of 
Eqs.  (A-1)  and  (A-8) . The  Intensity  autocorrelation  function  Rj^  for 
the  scatter  component  Is  given  In  terms  of  R^  and  from  Eq.  (A-16) , as 

Rj  (a"^)  - 1 - 4R^^(a"p)(i  - 4 ) + |Rg(Ap)|^  + Ib^Ca'p)!^  (A-31) 

S 'S' 

where 

- 1 - exp  Al|  (A-32) 

and 

R - i [R  + BR  ] . (A-33) 

XX  2 s s 
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FIGURE  A-1  CALCULATION  FOR  HIGHER-ORDER  PHASE  DERIVATIVES 
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When  Ap  ■0,  R (o)  - 1,  reduces  to  the  Intensity  scintillation 
2 

index  for  the  scatter  component. 


For  the  focus  component,  the  corresponding  autocorrelation  function 


is 


Rj.  (Ap)  = exp 


(A-34) 


where  <XfXf>  is  evaluated  by  using  Eq.  (A-30) . As  with  the  scatter  com- 
^ ^ 2 

ponent,  when  Ap  = 0,  = Ri^(o)  - 1. 

The  composite  intensity  autocorrelation  function  is  given  as: 


R (Ap)  = R^  R 

s ^f 


The  intensity  scintillation  index  is  given  as 


(A-35) 


(A-36) 


The  phase  scintillation  index  for  the  focus  component  is  given  approxi- 
2 2 

mately  by  o^AL  - o^.  This  approximation  is  valid  in  the  presence  of  weak 
focusing.  For  the  scatter  component,  a value  consistent  with  Eq.  (A-31) 
is  obtained  by  integrating  the  phase  distribution  function  as  described 
in  the  next  subsection. 


5.  Phase  Statistics  for  the  Scatter  Component — The  Tlatfield  Distribution 

If  x and  y are  jointly  Gaussian  but  otherwise  arbitrary  random  vari- 
ables, then 


1 


1 


f 

I 
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(A-37) 


/oa 

a P (a  C08  iji,  a atn  (^)  da 
xy 

u 


is  the  phase  probability  density  function,  where  P (...)  is  the  joint 

xy 

probability  density  for  the  quadrature  components  x and  y.  The  form  of 

P (...)  is  well  known,  and  it  has  been  shown  by  Hatfield  (Hatfield  and 
xy 

Rlno,  1975)  that  P^f^i)  admits  the  closed-form  solution 


2ii  a a ^1  - f?('J>) 

X V ^ 1 


(1  - p^)  exp 


2(1 


P ) 


+ fjC'^) 


V2n  - (1 


erfc 


■f2(<>) 


V2(l  - 

(A-38) 


where 


V5  2 2 2 

cos^  i()/p  - 2p  cos  4i  sin  (o  o ) + sin  iji/o  (A-39) 

X X y y 


^2^*^  - y^j  |cos  ip/a^  - p sin  <!)/ (o^Oy)| /f ^(41) 

+ y /sin  41/0^  - p cos  0/(o  o )/ /f , (4i) 
■'y  \ y X y J 1 


(A-40) 


and  a . 0 . and  0 are  the  individual  variances  of  x and  y,  and  0 is 
X y xy  ’^y 

their  covariance.  The  parameters  n and  n are  the  mean  values.  Finally, 

X y 


0/00 
xy  X y 


(A-41) 
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f 

I 

To  obtain  the  ms  phase  for  the  scatter  component,  the  Hatfield  dis- 
tribution, Eq.  (A-38) , Is  numerically  Integrated  in  Subroutine  CP  RMS. 

The  distribution  Itself  Is  evaluated  In  Subroutine  PHAS. 

6 . Two-Frequency  Correlation  Function 

The  two-frequency  correlation  function  for  Intensity  Is  a straight- 
forward extension  of  Eqs.  (A-31) , (A-3A) , and  (A-35)  In  Section  4 of  this 
appendix.  For  the  scatter  component  we  have 


s s 


4<x(^)x<2>> 

s s 


Af 


(A-42) 


i 


i 


! 


k 


where  the  superscripts  refer  to  the  frequencies  f^^^  - f - Af  and 
f^^^  “ f + Af.  The  correlation  functions  R (^;z)  and  B (^;2)  are  com- 

SI  SI 

puted  from  a generalization  of  Eq.  (A-17)  In  Section  2 of  this  appendix. 
It  can  be  shown  that 


(A-43) 


where  F(C)  is  computed  from  Eq.  (A-17b)  In  Section  2 with  Z replaced  by 


(A-44) 
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II 


In  Eq.  (A-44) , \ is  evaluated  at  the  center  frequency  f 
The  Input  correlation  functions  R and  B are  given  by  Eq. 
replaced  by 


(f(l)  f(2) 

(A-i6)  with  0 


)/2. 


Finally, 


(A-45) 


(A-46) 


For  the  focus  component,  Eq.  (A-30)  in  Section  3 is  used  with  X 
replaced  by  . The  correlation  coefficient  computed  Is 


'^1^2 


Ti 


<,(i)2y  4(2)2>*  j 


(A-47) 


which  varies  between  1 for  perfect  correlation  and  1^1  + 


* sf  >^)] 


-1/4 


[(1  - 


fffz 


^ yj  if  the  two  frequencies  are  completely  decorrelated . 

This  possibly  confusing  variation  occurs  because  the  average  value  of 
intensity  la  unity.  Independent  of  the  scintillation  level.  Thus,  C 

measures  the  correlation  between  intensity  fluctuations  at  the  two  fre- 
quencies rather  than  the  correlation  between  the  total  intensities 

7.  The  Empirical  RMS  Electron  Density  Model 

The  rms  electron  density,  AN,  is  calculated  in  RMSDN  as 


AN  - ANj,  4Nj^  + ANj^  + AN^ 
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(A-48) 


where  the  subscripts  refer  to  equatorial,  midlatitude,  high  latitude,  and 
auroral  respectively.  Each  term  Is  given  as  a function  of  mean  sunspot 
number,  R,  geomagnetic  latitude,  X^,  day  of  year,  D (1  ^ D £ 365),  time 
of  day  in  decimal  hours,  and  the  planetary  magnetic  index,  K^,  as  follows: 


AN,,  = 2.3  X 10  (1  + 0.04  R) 
£ 


' 

r 

exp  - 

(X^  - 10*)^ 
in 

+ exp 

(X  + 10“)^ 
m \ 

^ 2 

« 2 1 

(10*)^ 

(10*)'^  \ 

■ ^ 

/,  „ , r n(D  + 10)  . ir(D  + 10)D 

^1  - 0.4  [cos  ^ + yfr  cos  ^ 


exp 


(T  - 22.5) 
-2 


2 1 


el/m 


(A-49a) 


where 


- {1}  “ T ^ 22-5 


(A-49b) 


1.3  X 10 


*^1  + 0.33  cos  exp 


(X  - 32.5  )‘ 
m 


(15*)' 


el/m 


(A-50) 


ANjj  - 3 X 10^(l  + erf  el/m^ 


irT 


■) 


(A-51a) 


X._  - 68.5  - 1.8  K - (5  + 0.5  K ) cos  deg  (A-51b) 


- -p  p,  ^2 
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i 


■ 5 + 0.2  Xfa  deg 


(A-51c) 


3.2  X 10  K exp 
P 


— ^ el/m^  (A-51d) 

(1.1  Kp)2 


with 

R - Mean  'tunspot  number 

" Geomagnetic  latitude 

Xg  " Geographic  latitude 

D ” Day  out  of  365 

T “ Time  of  day  in  decimal  hours 

K • Planetary  magnetic  index. 

P 

All  parameters  in  the  model  refer  to  F- layer  penetration  point. 


In  the  subroutines,  MDLPRM,  additional  parameters  that  vary  with 
latitude  are  computed.  They  are  the  elongation  parameters  a and  b,  the 
outer  scale  for  the  focus  component,  a^,  the  layer  thickness,  L,  and  the 
layer  height,  H.  The  specific  formulas  are 


a 


(A-52) 


/ r\  - 20'n\ 

(l  . erf 


L - 200  -■  50  1 + erf 


km  (A-55) 


H - 500 


-4eerf[^]) 


km  . (A-56) 


The  remaining  scale  parameters,  and  a , are  given  In  terms  of  a,  as 

r 8 r 

- 0^/2. 5 (A-57a) 

a ■ a,/^  . (A-57b) 

s t 


8.  Coordinate  Transformations 

All  diffraction  calculations  are  performed  In  a downward-directed 
local  geomagnetic  coordinate  system.  The  receiver  coordinate  system  Is 
upward-directed  and  tangent  to  the  earth's  surface.  For  characterizing 
the  spatial  structure  of  the  field,  It  Is  necessary  to  convert  from  one 
coordinate  system  to  the  other.  The  two  coordinate  systems  are  shown  In 
Figure  A-2. 

The  coordinate  rotation  that  takes  the  receiver  systems  parallel  to 
xyz  Is 


(A-58) 
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where  x >•  (x,  y,  z)  , x = (x  , y , z ) , and  the  matrix  R has  elements 

P P P P 


R, , ” cos  A cos  + sin  A sin  A„ 
11  p R p R 


R,  _ “ cos  A sin  A„  + sin  A cos  A„)  cos 
12  p R p R R 


R,  _ = (cos  A sin  A - sin  A cos  A„)  cos  <p„ 
Xo  p K p K R 


21 


(-  sin  A cos  A„  + cos  A sin  A„)  sin  <p 
p R p R p 


(A-59a) 

(A-59b) 

(A-59c) 

(A-59d) 


'22 


[(sin  A sin  A„  + cos  A cos  A„)  sin  (p  sin  d„  + cos  <p  cos 

p R p R p R R 


(A-59e) 

* [-(sin  A sin  A + cos  A cos  A ) sin  <p  cos  <P  + cos  (p  sin  <0^] 
pR  pR  pR 

(A-59f) 
(A-59g) 
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(sin  A cos  A_  - sin  A_  cos  A ) cos  <p 
p R R p ^p 


R22  ■ [-(sin  A^  sin  A^j  + cos  A^  cos  A^)  cos  <p^  sin  + sin  <p^  cos  v’dI 
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P R' 

(A-59h) 

[(sin  Ap  sin  + cos  A^  cos  Aj^)  cos  <p^  cos  + sin  sin 

(A-59i) 


where  A denotes  latitude,  <p  denotes  longitude,  and  the  subscripts  p and 
R denote  penetration  point  and  receiver,  respectively. 

The  rotation  to  local  geomagnetic  coordinates  is 


X “ sin  6x  + cos  6y 
P P 

y ■ cos  6x  - sin  6y 
P P 


(A-60a) 

(A-60b) 
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* ' I 


-1  T 

where  6 Is  the  magnetic  declination  angle.  We  note  that  R » R . The 
superscript  T denotes  transpose.  To  characterize  the  spatial  autocor- 
relation function  at  the  receiver,  we  use  the  transformation 

Ax  = ^2’  ^*R  ^ ^^12  ^ ~ ^21  ^^R 

(A-61a) 

Ay  = cos  6 - sin  6)  Ax^^  + (Mj^2  ^ ~ ^22 

(A-61b) 

Az  = - AXj^  - M22  • (A-61c) 

The  family  of  ellipses  generated  by  f(Ap^)  = consi'.  Is  characterized 
by  the  orientation  angle 

= 1/2  arctan  ~J~cj  (A-62) 

measured  counterclockwise  from  the  x-axis  to  the  semi-major  axis  of  the 
ellipse.  The  axial  ratio  Is 

AR  = [(A  + C + D)/(A  + C - D)]*  (A-63) 

where 

D = [(A  - C)^  + B^]*  . (A-64) 

Since  we  are  using  a downward-directed  geomagnetic  coordinate  system,  the 

x-axis  points  in  the  magnetic  southward  direction.  To  calculate  the 

anisotropy  In  the  receiver  plane,  the  coefficients  C,  -B,  and  A In  Eqs. 

2 

(A-62),  (A-63),  (A-64)  are  replaced  by  the  coefficients  of  Ax  , Ax  Ay  , 

2 R R K 

and  Ay  , respectively. 

R 

9.  The  Nominal  Orbit  Code 

A satellite  In  circular  orbit  Is  easily  described  mathematically. 

The  relevant  parameters  are  shown  In  Figure  A-3.  Within  the  orbit  plane 
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PLANE  AND  EOUATORIAL  PLANE) 


LA-4259-21 


FIGURE  A-3  GEOMETRY  OF  SATELLITE  ORBIT 


the  satellite  position  is  determined  by  the  time,  t,  the  orbital  period, 
T,  and  a time  parameter,  t^,  which  determines  the  time  of  the  equator 
plane  crossing.  The  orbit  plane  itself  is  located  by  two  angles — namely, 
the  inclination  relative  to  the  north  celestial  pole,  9^,  and  the  loca- 
tion of  the  Orbit  plane  relative  to  the  sun,  h^,  which  we  take  to  be 
fixed. 

The  position  of  the  satellite  in  the  earth-centered  system  XYZ 
shown  in  Figure  A-3  is 

X - R (-  cos  X (t)  sin  9.  sin  9(t)  - sin  X (t)  cos  9(t))  (A-65a) 
s'  s i s 

Y - R (sin  X (t)  sin  9,  sin  9(t)  - cos  X (t)  cos  9(t))  (A-65b) 

s S ■ 1 = 


Z ■ R cos  9.  sin  9(t) 


(A-65c) 


i 


I 


i 


m 


where 

and 


A^(t) 


0(t)  - 


If  (t  - 12  - h^) 


^ (t  - t ) 
T O 


(A-66) 


(A-67) 


The  orbital  period  and  are  related  by  the  orbit  equation 


pi*  / 

T - 2Tr  / 

I m Gl  I 

L ® J / 
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(A-68) 


In  Eq . (A-68) , is  the  mass  of  the  earth  and  G is  the  gravitation 
constant. 

From  Eqs.  (A-65a)  through  (A-66),  one  can  show  that  the  latitude  and 
longitude  of  the  satellite  are  given  by  the  formulas 

2 2 I 

SLAT(t)  - arctan  |sln  8(t)  cos  0^,  [cos  0(t)  + sin  0^  sin  0(t)]} 

(A-69) 

SLON(t)  - arctan  I - (sin  0^  sin  6(t)  cos  X(t)  + cos  0(t)  sin  X(t))  , 
(cos  0(t)  cos  X(t)  - sin  0^  sin  0^(t)  sin  X(t))^.  (A-70) 

respectively.  In  subroutine  FINDORB,  Eqs.  (A-69)  and  (A-70)  are  solved 
Iteratively  for  the  orbital  parameters  t,  h^,  and  0^,  for  a given  period 
T and  a set  of  two  points  on  the  earth.  A reference  time  is  specified 
for  the  first  point. 


The  algorithm  allows  the  user  to  generate  a local  circular  orbit 
that  will  approximate  the  trajectory  of  any  actual  satellite,  at  least 
over  a small  Interval. 


Once  the  orbit  is  specified,  it  is  straightforward  to  calculate  the 
satellite  scan  velocity  v^,  which  is  usually  the  dominant  contribution 
to  the  time  structure  of  the  signal.  The  routine  is  useful  in  general 
for  generating  parameter  variations  that  are  representative  of  the 
various  typical  orbiting  satellite  conditions  a user  might  encounter. 
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